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The  MOS  capacitors  were  irradiated  with  2-Mev  a  particles  or 
with  700-kev  protons.  The  flatband  voltage  shift  was  determined 
by  using  the  well-known  capacitance-voltage  (C-V)  technique. 
This  voltage  shift  is  proportional  to  the  amount  of  ionization 
that  escapes  recombination  in  the  S: 

C 

The  experimental  results  agree  very  well  with  a  modified 
version  of  the  columnar  recombination  model  first  proposed  by 
George  Jaffd  in  1913.  Jaff€  developed  an  approximate  analytical 
solution  that  is  not  entirely  satisfactory.  However,  by  using  a 
large  digital  computer,  one  can  solve  Jaffa's  exact  differential 
-f  -niMatinn  numei'ically.jj  If  the  column  radius,  b,  is  taken  to  be 
3.5  nm,  the  result"^  of  a  finite  difference  calculation  fit  the 
experimental  results  very  well.  At  fields  of  2  MV/cm,  the  yield 
of  charge  is  about  10  percent  for  incident  a  particles  and  about 
25  percent  for  protons. 

These  results  are  much  less  than  the  yield  that  one  would 
expect  for  other  forms  of  radiation.  For  this  reason,  one  can 


probably  conclude  that  electronic  devices  will  probably  not  fail 
from  single  particle  induced  ionization  until  the  devices  are 


much  smaller  than  those  now  on  the  drawing  boards.  — 
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1 .  INTRODUCTION 


In  recent  years,  the  effect  of  a  particles  and  other  ions  passing 
through  semiconductors  has  been  a  subject  of  intense  interest.  In  1978, 
May  and  Woods showed  that  so-called  soft  errors  or  upsets  in  micro¬ 
electronic  circuits  were  sometimes  caused  by  a  particles  from  radio¬ 
active  contaminants  in  the  ceramic  packaging  around  the  circuit.  When 
such  a  particle  passes  through  a  silicon  (Si)  chip,  it  generates  a  large 
number  of  electron-hole  pairs.  Sometimes  enough  charge  is  collected  to 
upset  a  memory  element  (that  is,  to  change  its  state).  Since  this 
effect  was  identified,  several  other  investigators^"^  have  studied  the 
effects  of  other  kinds  of  cosmic  rays  and  other  charged  particles  in 
semiconductors . 

In  this  study,  we  consider  the  effects  of  heavy  charged  particles  in 
insulating  silicon  dioxide  (Si02)  films  rather  than  in  semiconductors. 
Amorphous  Si02  is  an  insulator  of  great  fundamental  interest  because  it 
has  a  high  dielectric  strength  and  its  electronic  properties  have  been 
studied  at  fields  an  order  of  magnitude  or  more  higher  than  most  other 
insulators.®  Thermally  grown  Si02  is  also  a  material  of  great  practical 
interest  because  of  its  role  in  metal-oxide-semiconductor  (MOS)  micro¬ 
electronic  circuits  and  in  integrated  optics.  It  has  long  been  known 
that  Si02  is  extremely  sensitive  to  ionizing  radiation,  and  its  response 
to  electron  beam  and  gamma  radiation  has  been  studied  extensively. 
References  7  to  1 1  are  a  representative  but  far  from  complete  bibli¬ 
ography.  Thus,  there  are  fundamental  reasons  as  well  as  practical 
reasons  to  study  the  effects  of  heavy  charged  particles  in  thin  films  of 
amorphous  SIO2. 

In  the  work  reported  here,  we  conducted  a  series  of  experiments  in 
which  MOS  capacitors  were  exposed  to  a  beam  of  ions,  and  the  fraction  of 
charge  escaping  recombination  was  determined.  Then  we  adapted  the  Jaff^ 
columnar  recombination  model  to  fit  a  curve  to  the  experimental  results 
with  reasonably  good  agreement.  The  model  can  be  used  to  calculate  the 
recombination  for  particles  and  energies  different  from  those  used  in 
the  experiments.  Both  the  experiments  and  the  model  are  discussed  in 
detail  in  this  report. 


2.  HISTORICAL  BACKGROUND 

The  problem  of  ionization  of  insulators  by  radiation  is  among  the 
oldest  in  modern  physics.  The  theory  of  ionization  and  recombination  in 
gases  was  first  based  on  the  assumption  of  uniform  ionization  throughout 
the  total  volume  of  the  gas.  Bragg  and  Kleeraan^^'^^  were  the  first  to 
show  that  for  a  particle  irradiation,  the  recombination  is  much  stronger 
than  can  be  explained  by  a  theory  based  on  uniform  ionization.  They 

*See  references  in  Literature  Cited  section. 
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both  tried  to  explain  their  results  in  terms  of  an  "initial  recombina¬ 
tion"  model  (also  now  called  geminate  recombination),  which  assumed  that 
ions  in  a  pair  originating  from  the  seune  molecule  have  a  strong  special 
tendency  to  recombine  with  each  other.  This  gemi..ate  recombination 
model  has  been  further  developed  by  Smoluchowskl , ^ Onsager,^^  and 
others, and  it  has  proved  useful  in  treating  many  problems, 
although  some  problems  seem  to  require  more  complicated  analysis. 
However,  Moulin^** showed  that  it  did  not  apply  to  a  particles  in 
gases . 

Langevln^^'^^  was  the  first  to  realize  that  the  ionization  produced 
by  a  particles  was  in  extremely  dense  columns  about  the  particle  track 
so  that  the  recombination  would  be  much  stronger  than  predicted  by  a 
theory  based  on  uniform  ionization.  Langevin  argued  that  the  recombina¬ 
tion  was  governed  by  the  following  equation: 


-cm^n 


(1) 


where  n^  is  the  density  of  positive  (or  negative)  charge.  He  also 
obtained  an  analytical  expression  for  the  recombination  coefficient  a. 
The  results  obtained  by  Moulin  were  consistent  with  the  Langevin 
hypothesis. 

The  Langevin  theory  was  refined  by  JaffS  in  a  classic  paper  in 
1913.^®  JaffS  added  a  diffusion  term  to  account  for  random  thermal 
motion  and  a  drift  term  to  account  for  the  effect  of  a  constant  external 
field  to  the  Langevin  equation  (1).  Thus,  the  evolution  of  the  columns 
of  charge  produced  by  an  a  particle  is  described  by  the  following 
differential  equation : 

3n  3n 

—  =  D^72n^  ^  '  <2) 


where  the  a  particle  travels  along  the  z-axis,  D  is  the  diffusion 
coefficient,  ii  is  the  mobility  of  the  charged  particles,  and  E  is  the 
component  of  the  field  in  the  x-dlrection.  D  and  p  may  be  different  for 
the  positive  and  negative  species.  Jaff£  was  not  able  to  solve  this 
equation  exactly,  so  he  neglected  the  recombination  term  and  solved  the 
rest  of  the  equation,  getting  the  standard  diffusion  result.  Then  he 
reintroduced  the  effect  of  the  recombination  term  by  letting  the  total 
number  of  charges,  N,  vary  with  time,  t.  By  this  procedure,  he  obtained 
the  following  approximate  expression  for  the  amount  of  charge  escaping 
recombination : 


l/2 


N«e 


n 


where 

y  is  the  yield  of  charge  (fraction  escaping  recombination), 

Mq  is  the  total  number  of  ion  pairs  per  centimeter, 
e  is  the  electron  charge, 

e  is  the  relative  dielectric  constant  of  the  medium, 

Cq  is  the  permittivity  of  free  space, 
b  is  the  effective  column  radius, 

E  is  the  total  applied  field, 

6  is  the  angle  between  the  field  and  the  particle  track. 

Since  Ng  and  b  were  not  known  in  1913,  Jaff4  treated  them  as  adjustable 
parameters,  and  he  was  able  to  obtain  good  fits  of  equation  (3)  to  his 
original  experiments  and  also  to  a  later  series  of  experiments. 

JaffS  did  most  of  his  experiments  with  gases,  although  he  also 
worked  briefly  with  organic  liquids.  Gerritsen^ did  a  series  of 
experiments  with  a  particles  passing  through  liquified  gases,  and  he 
obtained  results  that  did  not  agree  very  well  with  Jaffa's  approximate 
solution  (eq  3).  In  an  effort  to  explain  Gerritsen's  results, 
Kramers^ ^  developed  a  different  approximate  solution  to  Jaffa's 
differential  equation  (eq  2).  Kramers  neglected  the  diffusion  term  and 
obtained  an  analytical  solution  to  the  remaining  equation.  Then  he 
attempted  to  reintroduce  the  effect  of  the  diffusion  term.  This  pro¬ 
cedure  was  not  very  satisfying  since  the  diffusion  correction  is 
unphyslcal,  and  the  agreement  between  Kramers'  theory  and  Gerritsen's 
experiments  also  was  not  as  good  as  one  could  wish. 

One  can  see  from  the  preceding  discussion  that  the  ionization 
produced  by  a  particles  passing  through  matter  is  an  old  problem  that 
has  been  studied  extensively.  A  great  deal  has  been  learned,  but 
significant  questions  remain  even  for  materials  that  were  studied  in  the 
past.  Ftor  a  fundamentally  different  material  such  as  thermally  grown 
Si02,  very  little  is  known  with  certainty  until  experiments  are 
performed. 


3.  EXPERIMENTAL  PROCEDURE  AND  RESULTS 

In  our  experiments,  we  used  dry  oxide  MOS  capacitors  on  n-type  Si 
prepared  by  Hughes  Aircraft  Corp.  The  samples  had  an  oxide  thickness  of 
135  nm  and  an  electrode  area  of  0.0032  cm^  (0.025-ln.  dieuneter).  The 
electrodes  were  thin  aluminum  (105  nm  thick)  with  bonding  pads  on  the 
edge  for  electrical  contact.  The  samples  were  mounted  on  TO-5  headers, 
and  electrical  contact  was  made  by  thermocompression  gold  wire  bond. 

The  samples  were  irradiated  under  positive  bias  at  approximately 
liquid  nitrogen  temperature  (~77  K).  The  samples  were  mounted  in  a 
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copper  (Cu)  cold  finger,  which  was  In  contact  with  a  liquid  nitrogen 
reservoir.  The  temperature  was  monitored  with  a  Cu-constantan  thermo¬ 
couple  , 

When  an  ion  or  any  Ionizing  radiation  passes  through  the  SIO2  film, 
it  loses  energy  by  creating  electron-hole  pairs.  At  liquid  nitrogen 
temperature,  the  electrons  are  swept  out  of  the  Si02  almost  instan¬ 
taneously,  but  the  holes  have  only  a  negligible  mobility.  For  excunple, 
Hughes  has  shown^^  that  the  electrons  have  a  mobility  of  20  cm^/V-s  at 
room  temperature  and  40  cm^/V-s  extrapolated  to  liquid  nitrogen  temper¬ 
ature.  Thus,  for  fields  on  the  order  of  10^  V/cm  and  an  SIO2  thickness 
on  the  order  of  10“^  cm,  all  the  electrons  are  swept  out  of  the  Si02  in 
less  than  1  ps  both  at  room  temperature  and  at  77  K. 

The  holes,  on  the  other  hand,  have  a  much  lower  mobility,  and  they 
do  not  move  by  a  conventional  "drift"  process.  The  holes  move  by  a 
dispersive  hopping  process,  which  has  been  extensively  studied^' ' 35 
and  which  has  been  shown  to  be  both  temperature  dependent 3^  and  field 
dependent. 3^  At  liquid  nitrogen  temperature,  the  mobility  of  the  holes 
is  10~^^  cm^/V-s  or  less,  depending  on  the  applied  field.  It  has  been 
shown®  that  hole  transport  is  negligible  at  times  of  100  s  or  even  more 
at  applied  fields  of  2  MV/cm  or  less.  At  an  applied  field  of  3  MV/cm, 
the  transport  is  not  a  large  effect,  but  it  is  detectable.  Thus,  the 
electrons  are  removed  from  the  Si02  "instantaneously,"  and  the  holes  are 
frozen  in  place  for  times  long  compared  with  the  experimental  times.  At 
room  temperature,  on  the  other  hand,  significant  hole  transport  would 
occur  during  the  experiment,  and  the  results  would  be  difficult  to 
interpret . 

The  a  particle  and  proton  irradiations  were  carried  out  at  the  Naval 
Surface  Weapons  Center  (NSWC)  Van  de  Graaff  facility.  The  a  particles 
had  a  kinetic  energy  of  2  MeV,  and  the  protons  had  a  kinetic  energy  of 
700  keV.  In  both  irradiations,  the  beam  was  brought  into  the  test  cham¬ 
ber  and  scattered  by  a  200-nm  Cu  foil.  The  particle  flux  was  determined 
by  using  two  surface  barrier  detectors :  one  was  in  the  sample  position 
(about  6  deg  above  the  beam  axis),  and  the  second  (monitor)  detector  was 
off  at  a  large  angle  (>45  deg)  from  the  beam  axis  (fig.  1).  The  ratio 
of  counts  in  the  two  detectors  was  nearly  constant  in  a  series  of  cali¬ 
bration  exposures  (table  1). 

In  principle,  the  ratio  of  counts  in  the  two  detectors  could  be  cal¬ 
culated  from  the  Rutherford  cross  section  if  the  angle  were  measured 
precisely.  But  if  the  seunple  were  raised  0.05  in.  (1.27  mm),  that  is, 
if  the  angle  were  Increased  less  than  0.5  deg,  the  ratio  of  counts  in 
the  two  detectors  would  be  reduced  approximately  30  percent.  Since  our 
test  capacitors  were  located  anywhere  on  a  TO-5  header,  that  is,  not 
necessarily  centered,  we  could  get  large  variations  in  particle  count  at 
the  s^unple  for  a  constant  particle  count  at  the  monitor  detector.  Few 
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Figure  1 .  Experimental  schematic  for  o  par¬ 
ticle  and  proton  exposures. 


TABLE  1 .  RESULTS  OF  CALIBRATION  EXPOSURES  OF  a  PARTICLES 
AND  PROTONS 


Calibration 

Exposure 

Sample 

position 

Monitor 

detector 

Exposure 

(s) 

Initial^ 

1 

911769 

11032 

138.2 

hallum  lona« 

2 

856319 

10023 

90.4 

l-inm  apar* 

3 

858821 

9985 

48.5 

ture 

4 

850685 

10017 

50.2 

5 

858404 

10022 

80.2 

Sample 

1 

607824 

9995 

51.6 

raised 

2 

606287 

9988 

53.7 

1 *27  mm 

3 

598782 

10028 

51.9 

4 

594426 

10027 

51.2 

5 

593954 

10025 

52.8 

Proton, 

t 

501689 

10063 

1-imn 

2 

503328 

10050 

• 

aperture 

3 

497349 

10051 

. 

4 

503061 

10073 

. 

5 

510102 

10074 

- 

things  in  nuclear  physics  are  as  well  established  as  the  Rutherford 
cross  section,  but  it  varies  so  strongly  with  position  at  small  angles 
that  we  could  not  rely  on  it  to  determine  the  incident  particle  count  at 
the  sample. 

The  method  finally  settled  on  for  counting  Incident  particles  was  to 
use  the  sample  as  its  own  detector*  After  the  sample  had  been  exposed 
at  positive  bias  and  the  flatband  voltage  shift  (AVp^)  had  been  deter¬ 
mined,  the  bias  was  reversed  so  that  the  sample  was  in  depletion.  The 
NOS  capacitor  was  then  used  like  a  surface  barrier  detector  with  the 
depletion  region  in  the  Si  being  the  sensitive  volume.  Then  the  expo¬ 
sure  was  repeated,  and  the  current  pulses  in  the  capacitor  were  counted 
by  a  multichannel  spectirum  analyzer* 


Figure  2  is  a  schematic  diagreun  of  a  sample  during  the  particle 
counting.  As  an  a  particle  passes  through  an  MOS  capacitor,  it  creates' 
electron-hole  pairs  in  the  Si02  and  in  the  Si  substrate.  The  dashed 
line  represents  the  depth  in  the  Si  of  the  depletion  layer.  The  free 
charge  in  the  depletion  layer  appears  as  a  current  pulse,  which  is 
counted  by  the  spectrum  analyzer.  A  typical  spectrum  is  shown  in  figure 
3.  A  background  measurement  Indicated  that  the  low-energy  counts  (below 
the  marker)  were  present  with  the  beam  blocked  off,  but  only  about  1000 
counts  above  the  mark  were  noise. 

The  aluminum  gate  electrodes  are  thin  enough  that  the  energy  lost  by 
the  charged  particles  is  almost  negligible.  A  2-MeV  a  particle  passing 
through  105  nm  of  aluminum  (Al)  loses  slightly  less  than  0.04  MeV.^® 
For  the  protons,  the  energy  lost  in  the  gate  electrode  is  even  less, 
0.005  MeV.  Thus,  a  known  number  of  particles  of  known  energy  were  inci¬ 
dent  on  the  Si02  film.  The  amount  of  charge  created  could  then  be  cal¬ 
culated  very  accurately.  The  amount  of  charge  that  escaped  recombina¬ 
tion  was  determined  by  measuring  the  AVp3  using  the  well-known 
capacitance- voltage  (C-V)  technique.®®  Typical  C-V  curves  are  shown  in 
figure  4.  The  irradiation  took  place  at  77  K,  and  the  preirradiation 
and  postirradiation  C-V  traces  were  recorded.  The  average  fluence  of  a 
particles  and  protons  on  the  sample  was  1.6  x  10®/cm^  and  1.06  x 
10®/cm^,  respectively.  Fbr  the  very  thin  Si02  (135  nm) ,  the  energy 
loss  per  unit  path  length  was  uniform  across  the  Si02  for  these  high 
energy  particles.  The  calculated  total  dose  in  the  Si02  was  3.8 
krad  (Si02)  from  the  a  particles  and  0.54  krad  (Si02)  from  the  protons. 


a 


Figure  2.  Sample  schematic  showing  metal,  Si02,  and 
depletion  region  of  Si  substrate. 


LOG  OF  PARTICLE  COUNT 


Figure  3.  Typical  a  particle  spectrum  as  recorded  by  metal- 
oxide-semiconductor  capacitor;  logarithm  of  parti¬ 
cle  count  is  plotted  against  particle  energy;  low- 
energy  noise  (below  mark)  is  subtracted  out. 


Figure  4.  Typical  capacitance-voltage  curve 
from  this  experiment. 
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The  total  AVpg  was  assumed  to  be  proportional  to  the  average  yield 
of  holes  that  escaped  recombination  and  were  immobilized  in  the  Si02* 
The  AVp3  is  proportional  to  the  first  moment  of  the  charge  distribu¬ 
tions,  and  in  this  case  the  holes  are  created  and  frozen  in  place  uni¬ 
formly  across  the  Si02>  tinder  these  conditions. 


AV 


FB 


1 

C 

OX 


(4) 


reduces  to 


AV 


FB 


1-2- 


2  C 


ox 


(5) 


where 

Cqjj  =  oxide  capacitance,  eeo/iox» 
e  =  3.85  for  Si02* 
eo  =  8.85  *  10"^**  f/cm, 
fox  =  oxide  thickness, 

Q  =>  charge,  f(E)Qo, 

f(E)  =  experimentally  determined  fractional  yield  of  charge  sur¬ 
viving  recombination, 

Qo  =  total  amount  of  positive  charge  generated  by  irradiation. 


fox  I’S  *  10”^^  coulomb/hole 
cos  6  18  eV 


(6) 


where 

dT/dx  “  energy  loss  per  unit  path  length, 

18  eV  -  electron-hole  pair  creation  energy  that  has  been 
reported.**®”**^ 


14 


.  f 


The  procedure  was  to  measure  the  AVp3  as  a  function  of  the  applied 
field.  The  fractional  yield  was  then  determined  by  dividing  the  meas¬ 
ured  AVpb(E)  by  the  calculated  maximum  AVpg,  assuming  no  recombination. 

Since  one  would  expect  the  charge  yield  to  depend  on  the  angle 
between  the  field  and  the  particle  beam,  we  performed  three  series  of 
measurements:  (1)  with  a  particles  Incident  at  45  deg;  (2)  with  a 
particles  as  close  as  possible  to  normal  Incidence,  6  deg;  and  (3)  with 
protons  Incident  at  45  deg  (tables  2  to  4;  fig.  5  and  6). 


TABLE  2.  HELIUM  lOM  RESULTS  AT  45  DEG 


Shot 

Sample 

Applied 
field 
(MV/ cm) 

Change  in 
threshold 
voltage 
(V) 

a 

particle 

count 

Monitor 

count 

Fractional 

yield 

1 

118 

1,46 

0,106 

499062 

25100 

0.087 

2 

159 

1.11 

0.069 

423646 

25534 

0.066 

3 

104 

0,89 

0,051 

363908 

25025 

0.058 

4 

153 

0.37 

0.031 

355028 

24993 

0.036 

5 

125 

0.22 

0.026 

339447 

25108 

0.031 

6 

160 

0.15 

0.028 

469308 

25041 

0.024 

7 

137 

0.07 

0.015 

393174 

25000 

0.016 

S 

129 

0 

0.020 

410312 

25041 

0.014 

9 

112 

0,37 

0.033 

370939 

25081 

0.036 

10 

30 

0.52 

0,045 

349000 

24822 

0.055 

11 

130 

0.74 

0,047 

349000 

24822 

0.055 

12 

152 

2.22 

0.055 

182605 

25049 

0.113 

13 

135 

1.85 

0,098 

418907 

25052 

0,096 

TABLE  3. 

HELIUM  ION 

RESULTS 

AT  6  DEG 

Shot 

Sample 

Applied 

field 

(MV/cm) 

Change  in 
threshold 
voltage 
(V) 

O 

particle 

count 

Monitor 

count 

Fractional 

yield 

1 

132 

1.48 

0.083 

735248 

25074 

0.065 

2 

146 

0.74 

0.051 

684609 

25055 

0.041 

3 

113 

0.37 

0.024 

589315 

25034 

0.024 

4 

161 

1.11 

0.057 

499087 

25062 

0.050 

5 

102 

1.85 

0,094 

498744 

25126 

0.077 

6 

96 

2.22 

0.102 

644746 

25083 

0.091 

7 

126 

2,59 

0.071 

520040 

25072 

0.095 

8 

101 

2.96 

0.114 

666302 

25084 

0.098 

9 

140 

0.15 

0.029 

574529 

25033 

0.020 

Applied 

field 

(MV/cn) 


Chancja  in 
threshold 
vol tage 
(V) 


Particle 

count 


Monitor 

count 


fractional 

yield 


electric  FIELD  (MV/eill) 

Figure  6.  Experimental  results  for  protons  inci¬ 
dent  at  45  deg* 


One  can  see  from  figures 
5  and  6  that  the  fraction  of 
charge  escaping  recombina¬ 
tion  is  only  a  small  frac¬ 
tion  of  the  total  charge 
generated  by  a  highly  ion¬ 
izing  particle.  Par  com¬ 
parison,  we  show  in  figure  7 
the  results  of  these  experi¬ 
ments  along  with  previously 
published  results  for  elec*- 
trons.^^'**®  Fbr  high-energy 
electrons,  the  yield  is 
close  to  100  percent  at 
fields  above  ~2  MV/cm;  but 
for  a  particles  the  yield  is 
on  the  order  of  10  percent 
at  2  MV/cm,  and  for  protons 
the  yield  is  ~25  percent  at 
2  MV/cm.  The  results  in 
figure  7  show  dramatically 
the  effect  of  recombination 
along  the  track  of  a  heavy 
ion. 


ELECTRIC  FIELD  (MV/cm) 


Figure  7.  Comparison  of  a  particle  and 
proton  results  with  previously 
published  results  for  elec¬ 
trons  . 
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4.  ANALYSIS 


This  section  discusses  the  physical  model  that  describes  the  experi¬ 
mental  results*  In  particular,  we  detail  the  Jaff4  model^^  as  orig¬ 
inally  developed  and  as  applied  here* 

JaffS  began  by  assuming  that  the  charge  distribution  produced  by  an 
a  particle  is  uniform  along  the  z-axis,  that  is,  along  the  particle 
track*  The  initial  density  of  free  charge  about  the  axis  is  assumed  to 
be  Gaussian  and  given  by 


n^(r,o) 


**0  -rVb2 
-  e 

irb^ 


(7) 


where  Nq  is  the  total  number  of  electron-hole  pairs  produced  per  unit 
path  length  and  b  is  the  column  radius* 

Then  Jaff4  treated  the  case  in  which  no  external  field  is  present 
and  the  evolution  of  the  charge  distributions  is  given  by 

3n 

=  DV^n^  -  om^n_  ,  (8) 


where  the  first  term  on  the  right  is  the  diffusion  term  and  the  second 
term  is  the  recombination  term*  Jaff4  chose  to  develop  an  approximate 
solution  to  equation  (8)  by  neglecting  the  recombination  term  and 
solving  the  diffusion  equation  first*  He  obtained  the  standard  result 

.  ”o  -r2/f4Dt+b21 
n  (r,t)  =  - -  e  '  ■'  *  (9) 

*  n(4Dt  +  b2) 

Then  he  allowed  for  the  effect  of  recombination  by  replacing  Nq  with 
N(t)  in  equation  (9)  and  substituting  equation  (9)  into  equation  (8)* 
The  result  was  an  expression  for  dN(t)/dt,  which  could  be  Integrated  to 
obtain  the  result 


or 


N(t) 


.  4Dt  +  b^ 


b2 


n. (n,t) 

X 


^-r2/(4Dt+b2} 


.  .  ,  4Dt  +  b2  irf4Dt  +  b21 

1  +  -T-T  In  - -  '■ 


8TrD 
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(10) 


(11) 


.  r 


Next,  Jaffe  calculated  the  fraction  of  charge  diffusing  outside  a 
radius ,  R  ( comparable  to  the  separation  between  columns ) ,  which  was 
considered  to  have  escaped  recombination  at  any  time,  t.  This  calcu¬ 
lation  was  performed  by  integrating  equation  (11)  over  space  and 
with  the  fraction  of  charge  surviving  at  infinite  time, 

where 


4Dt  +  b^ 


°  b2 


dc  =  dt  . 

(4Dt  +  b2)2 


Since  R  is  generally  large  compared  with  b,  it  is  sufficient  to 
evaluate  the  integrand  at  the  upper  limit,  t  The  integrand  in 

equation  (12)  is  plotted  in  figure  8,  and  the  integral  can  be  evaluated 
numerically.  The  yield  for  2-MeV  a  particles  in  Si02  is  less  than  1 
part  in  10^. 


The  zero  field  case  corresponds  to  two  cylinders  of  charge,  one 
positive  and  one  negative,  spreading  out  under  the  influence  of  diffu¬ 
sion  and  undergoing  recombination  at  the  same  time.  Next,  Jaff€ 
considered  these  same  cylinders  in  the  presence  of  an  external  electric 
field.  The  basic  assumptions  here  are  that  (1)  the  charge  densities  are 
sufficiently  small  that  the  electron-hole  interactions  are  small  com¬ 
pared  with  interactions  with  the  field,  that  is,  e2/<r^3>  <<  eEX,  and 
(2)  screening  is  negligible  (b  <  ^pebye)*  first  of  those  conditions 
is  not  really  satisfied  until  considerable  recombination  has  taken 
place,  so  one  would  expect  the  theory  to  make  reasonable  predictions 
only  at  late  time. 

First,  Jaff4  considered  the  case  where  the  field  is  parallel  to  the 
column  axis  and  the  cylinders  move  as  indicated  schematically  in  figure 
9.  If  field  E  is  applied  between  the  plates  of  a  capacitor  with  sepa¬ 
ration  L,  then  the  time  for  the  cylinders  to  move  past  each  other  is 
given  by  T  *  L/2|jE.  (2m  is  replaced  by  +  M_  if  the  positive  and 
negative  charges  have  different  mobilities.)  The  region  where  the 
cylinders  overlap  is  described  by  the  zero  field  calculation  stated 


L06dN(l)/llt 


above,  but  no  further  recombination  takes  place  in  the  parts  of  the 
cylinders  where  the  positive  and  negative  distributions  no  longer 
overlap.  The  number  of  charges  of  each  sign  that  escape  from  the  volume 
of  overlap  in  time  increment  dt  is  2uEtl(t)  dt,  and  the  total  number  of 
free  charges  that  escape  recombination  is  given  by 


=  2viEN 


r _ 

“Jo 

8irD 


dt 


In 


4Dt  +  b^ 
b2 


(13) 


Using  the  logarithmic  integral. 


li(x) 


tabulated  by  Jahnke  and  Kmde,****  one  can  compute  the  fractional  yield, 
Y  =  N2/NqL,  which  can  be  shown  to  be 


Y 


^  y^e’^1(li  e^2  -  li  e^^) 


(14) 


where 


1 


SifD  .  ,  4DT  +  b^ 
y,  =  +  In  - 

^  “**0  b2 


Evaluating  equation  (14)  for  Si02  leads  to  the  results  pl'>tted  in  figure 
10  for  two  values  of  L. 

Second,  Jaff4  considered  the  case  where  the  field  is  perpendicular 
to  the  particle  track  as  indicated  schematically  in  figure  11.  He 
solved  the  following  differential  equation  by  first  neglecting  the 
recombination  term: 


3n 

?  y^E  3^  -  on^n^ 


(15) 
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Figure  10. 

Yield  calculated  for  nor¬ 
mally  incident  o  particle 
using  Jaffa's  procedure 
for  Si02. 


The  particle  track  is  here  taken  to  be  along  the  z-axiSf  and  the  field 
is  taken  to  be  along  the  x-axis.  When  this  equation  is  solved  without 
the  recombination  term,  one  obtains  the  result 


n 


± 


N 

'?r(4Dt  +  b^) 


exp 


(x  +  uEt)^  + 
4Dt  +  b2 


(16) 


where  r^  has  been  replaced  with  x^  +  y^.  Equation  (16)  is  the  standard 
diffusion  result  for  two  cylinders  of  charge  spreading  through  diffusion 
and  being  pulled  apart  as  shown  in  figure  12.  Then  Jaff4  substituted 
equation  (16)  into  equation  (15),  letting  N  be  a  function  of  t  (only), 
and  obtained  the  result 


N(t') 


r  exp(- 

I  \  4Dt  +  b^y 

•'0  4Dt  +  b^ 


1 

^  2ir 


dt 


(17) 
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ELECTRIC  FIELD 


Figure  11.  Schematic  of  two  cylinders  of  charge  moving  under 
influence  of  normal  field  and  parallel  field 
(that  is,  arbitrary  angle  of  incidence). 


ELECTRIC 

FIELD 


Figure  12.  Schematic  of  two  cylinders  of  charge  moving  under 
influence  of  normal  field. 
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If  equation  (17)  is  substituted  into  equation  (16),  one  obtains  the 
approximate  solution  to  equation  (15).  To  get  the  yield  of  charge,  one 
has  to  perform  the  integration  in  equation  (17)  from  t  =  0  to  t  =  “>, 
with  the  following  result : 


where 


S(z)  -  , 


2D^ 


(18) 


(19) 


Here,  ^  is  the  Hankel  function  of  the  first  kind. 

Third,  Jaff6  treats  the  case  of  arbitrary  angle  between  field  and 
particle  based  on  the  two  special  cases  already  considered  (fig.  13). 
If  the  angle  between  the  field  and  the  particle  track  is  taken  to  be  6, 
then  the  time  for  the  cylinders  to  move  past  each  other  is 

T'  =  L/2yE  cos  0  , 

which  is  analogous  to  the  parallel  field  case.  The  number  of  charges 
escaping  recombination  is 


r T' 

Nj  =  2mE  cos  0  N(t)  dt 


as  before,  where  E  has  been  replaced  by  the  relevant  component, 
E  cos  0.  The  fractional  yield  is  given  by 


y(E) 


(20) 
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particles  at  45  deg  compared 
with  Jaff4*s  approximate  solu¬ 
tion,  where  column  radius,  b, 
is  taken  to  be  adji.'atible . 

In  equation  (20),  one  can  recognize  equati'*?*.  (17)  v*th  field  E  replaced 
by  the  normal  component,  E  sin  0.  Equation  (20)  can  be  evaluated  with 
the  result 


y(E) 


1  + 


°^0 

8nD 


S(z*  ) 


(21) 


where  z*  *  sin^  0)/2D^.  Equation  (21)  is  only  approximately 
equal  to  equation  (20)  because  an  end  contribution  has  been  neglected, 
but  Jaff€  argues  that  it  is  small  in  practice.  For  fields  of  practical 
interest  in  Si02  (that  is,  high  fields),  one  can  use  the  asymptotic 
expansion  iH^^’(iz)  »  e"*/(TTx/2)l'2  and  get 

^<E)  “  |l  +  ^  [(|)^^\wE^sin  ejj  • 
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The  recombination  coefficient  a  was  determined  by  Langevin^®  to  be 

a  =  ye/cEfl  , 


and  equation  (22)  becomes 


Y(E) 


■  t  •  (!)''■ 


1-1 


4neeQbE  sin  6 


(23) 


which  is  identical  with  equation  (3).  Jaff4  was  able  to  explain  his 
experiments  adequately  with  this  model/  treating  both  Nq  and  b  as 
adjustable  parameters. However,  Nq  can  be  obtained  Independently 
by  using  the  Bethe  theory**^'**^  or  (recently)  Ziegler's  compilation,^®  so 
a  modern  researcher  has  only  one  free  parameter,  the  column  radius,  in 
applying  this  model.  In  figure  13,  we  replot  our  experimental  results 
for  a  particles  along  with  the  Jaff4  solution  for  several  values  of  b. 
One  can  see  that  the  agreement  between  the  JaffS  theory  and  the 
experiment  is  qualitative  at 
best.  The  curves  do  not  have 
the  same  shapes,  but  for  a 
narrow  range  of  E  values  one 
could  pick  a  value  of  b  that 
would  give  approximately  the 
correct  yield. 

In  measuring  the  yield  of 
charge  in  liquified  gases,  Ger- 
ritsen  encountered  similar 
problems . ®  ^ ^ 2  That  is,  the 
Jaff€  model  fitted  his  results 
only  roughly.  Kramers,®®  a  col¬ 
league  of  Gerrltsen,  proposed  a 
modification  of  the  Jaff€  treat¬ 
ment  to  try  to  explain  Gerrit- 
sen's  results.  Kramers  first 
neglected  the  diffusion  term  in 
equation  (15)  ( or  eq  2 )  and  then 
by  a  complicated  procedure  ob¬ 
tained  an  analytical  solution  to 
the  remaining  equation.  The 
result  is  plotted  in  figure  14 
along  with  the  Jaff4  solu¬ 
tion  for  comparison.  In  fig-  moucEO  fkld. i 

ure  14,  the  fractional  yield 

is  plotted  as  a  function  Figure  14.  Kramers'  approximate  solu- 
of  reduced  field,  f  =  P/Fq,  tion  compared  with  Jaffa's 

where  F  is  the  applied  field  and  expression. 
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Fq  is  a  reference  field.  When  Kramers  reintroduced  the  effect  of  the 
diffusion  term,  he  was  left  with  an  arbitrary  integration  constant,  c, 
which  has  no  physical  meaning.  The  dashed  lines  in  figure  14  are  the 
fractional  yield  for  two  selected  values  of  this  constant  C.  One  can 
see  from  figure  14  that  the  yield  is  rather  sensitive  to  the  value 
chosen  for  this  parameter. 


This  situation  is  extremely  unsatisfying.  It  is  one  thing  to  esti¬ 
mate  a  physical  quantity  such  as  the  column  radius,  which  is  unknown  but 
which  can,  in  principle,  be  determined  independently.  It  is  another 
thing  to  pluck  out  of  thin  air  a  constant  that  has  no  meaning  and  that 
cannot  be  checked,  but  that  nevertheless  is  critical.  Thus,  we  do  not 
use  Kreuners'  result.  But  we  should  note  in  passing  that  even  with  the 
diffusion  "fudge  factor,"  the  Kramers*  theory  did  not  fit  Gerrltsen's 
experiments  well.  In  the  end,  Gerrltsen  had  to  assume  that  the  column 
radius  took  on  different  values  at  high  fields  and  at  low  fields  to  get 
reasonable  agreement.  This  assumption  is  best  greeted  skeptically. 


The  value  of  Kramers'  discussion  is  that  he  points  out  that,  at 
least  at  low  temperatures,  the  recombination  term  in  equation  (15)  (or 
eq  2)  is  by  far  the  largest  term,  with  the  drift  term  next  and  the 
diffusion  term  smallest.  But  rather  than  neglecting  a  term  to  simplify 
the  analysis,  one  can  use  a  large  digital  computer  to  solve  the  exact 
equation  numerically.  Today,  even  a  pedestrian  researcher  can  perform 
analysis  unimaginable  to  either  Jaff4  or  Kramers  before  the  age  of 
microelectronics . 


The  approach  here  is  to  use  the  initial  Gaussian  distribution  postu¬ 
lated  by  Jaffi  (eq  7)  and  to  calculate  the  distribution  of  charges  using 
a  two-dimensional  finite  difference  code  and  equation  (15).  The  results 
can  be  generalized  to  three  dimensions  by  a  process  similar  to  that  used 
to  get  from  equation  (18)  to  equation  (21). 

The  code  written  for  this  analysis  calculates  the  initial  density  of 
charge  at  each  point  in  a  two-dimensional  Cartesian  grid  with  a  maximum 
of  200  x-values  and  100  y-values.  The  particle  is  assumed  to  travel 
along  the  z-axls,  and  the  field  is  assumed  to  be  along  the  x-axls. 
(Since  the  problem  is  symmetric  in  y,  it  serves  no  purpose  to  grind  out 
the  calculation  in  the  lower  half  plane.)  The  grid  spacing  la  variable 
in  principle  but,  for  all  the  calculations  reported  here,  the  points 
were  0.5  nm  apart.  Since  the  calculation  is  two  dimensional,  it  de¬ 
scribes  only  the  region  of  overlap  in  figure  11,  and  only  the  normal 
component  of  the  field  enters  the  problem.  Then  a  procedure  similar  to 
Jaffa's  is  necessary  to  generalize  the  result  to  arbitrary  angles.  That 
is,  one  can  follow  the  procedure  leading  to  equation  (13),  except  that 
the  Integral  is  evaluated  numerically. 
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The  code  calculates  the  positive  and  negative  charge  densities  from 
equation  (7)  for  each  point  in  the  grid  at  t  =  zero*  Then  it  calculates 


j) 


+  l»j)  +  n^{i  -  1,j)  +  n^(i,j  +  1) 
- 


n^(l,j  -  1)  -  4n^(i,j) 
h2 


(24) 


where  i  is  the  x  index,  j  is  the  y  index,  and  h  is  the  grid  spacing. 
Then  the  code  calculates  3n/3x  at  each  point  from  the  relation 


3n^(i, j) 
3x 


+  1»j)  -  n^(i  -  lo) 

2h 


(25) 


The  derivatives  3n^/3t  are  calculated  at  each  point  by  using  equation 
(15).  Then  the  charge  densities  are  updated  according  to  the  relation 

3n  (i, j,t) 

n^(i»j»t  +  At)  =  n^(i,j,t)  +  — -  At  ,  (26) 

where  At  has  to  be  chosen  with  some  care  to  keep  the  problem  from 
becoming  unstable.  After  the  new  densities  have  been  calculated,  an 
integration  routine  adds  up  the  surviving  charge,  and  N(t)  is  printed 
out.  In  addition,  Jg  N(t')  dt’  is  calculated  after  each  cycle  so  that 
the  problem  can  be  generalized  to  arbitrary  angles  of  incidence.  Then 
the  process  is  repeated  by  using  the  new  densities  and  usually  a  recal¬ 
culated  time  step. 


This  code  is  extremely  simple  in  that  it  consists  of  perhaps  150 
FORTRAN  statements.  However,  it  requires  large  amounts  of  memory  and 
processing  time.  The  grid  200  by  100  means  20,000  grid  locations  and 
there  are  eight  variables  (n^,  n_,  V^n^,  v2n_,  3n^/3x,  3n_/3x,  3n  /3t, 
and  3n_/3t)  calculated  at  each  point.  Thus,  there  are  ~160,000 
variables  calculated  per  cycle  for  several  hundred  cycles.  On  an  IBM 
System/370  Model  168  computer,  this  code  requires  about  1.4  Mbytes  and 
perhaps  10  to  15  min  of  run  time.  (For  a  complete  source  listing,  see 
appendix  A. ) 
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In  the  calculations  presented  here,  the  following  constants  were 
used : 


Nq  =  1.33  X  10®  cm“^  for  a  particles  and  2.97  x  10^  cm"^  for 
protons , 

=  10“®  cm^/V-s, 

=  40  cm^/V-s, 

a  =  (u+  +  U-)e/e£o  “  1*88  x  io~®  cm®/s, 

D+  =  p+kT/e  =6.5  X  lo~®  cm^/s, 

D_  =  p_kT/e  =  0.261  cm^/s, 

b  =  3.5  nm,  free  parameter  chosen  to  produce  agreement  with 
experiment. 


The  distribution  of  charges  at  t  =  0  is  shown  in  fiMre  15.  In  figures 
16  to  21,  the  distributions  are  shown  for  t  =  10“^®,  10"^**,  3  x  10”^**, 
10“^®,  and  3  x  10“^®  s  for  a  perpendicular  field  of  10®  V/cm.  One  can 
see  that  the  positive  charges  do  not  move,  and  most  of  the  charge  recom¬ 
bines  before  the  field  separates  the  negative  charges  from  them.  For  a 
much  smaller  field,  E  =  10**  V/cm,  the  distribution  at  t  =  3  x  10“^®  s  is 
shown  in  figure  22.  One  can  see  that  diffusion  of  the  negative  charges 
is  important  because  the  distribution  has  spread  out.  The  field  has 
moved  the  negative  charges  only  slightly  to  the  left,  but  the  recom¬ 
bination  has  eaten  away  both  charge  distributions  in  the  region  where 
they  still  overlap.  Figures  23  to  30  show  the  charge  distributions 
after  they  have  been  separated  or  nearly  separated  at  different  fields 
for  incident  a  particles.  One  can  see  from  the  figures  how  the  yield 
decreases  with  decreasing  field. 


These  results  are  summarized  in  figure  31,  in  which  the  solid  line 
joins  the  calculated  yield  values  for  oi  particles  incident  at  45  deg  and 
the  experimental  points  are  shown  for  comparison.  One  can  see  that  the 
agreement  between  theory  and  experiment  is  very  good  if  one  assumes  a 
column  radius  of  3.5  nm,  especially  at  high  fields.  At  low  fields,  the 
measured  yield  is  slightly  higher  than  calculated,  however. 
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This  difference  at  low  fields  is  probably  due  to  the  form  of  the 
diffusion  coefficient  that  was  used.  D  was  calculated  from  the  relation 


where  p  was  measured.  This  relation  is  valid  only  for  relatively  dilute 
solutions.  In  this  case/  the  maximum  ionization  density  at  t  =  0  is  ex¬ 


tremely  high,  between  10^“  and  10**  charges/cm^ 


Thus,  some  form  of 


concentration  dependent  diffusion  model**  ^  or  a  hot  electron  model  prob¬ 
ably  is  appropriate.  Such  a  model  would  increase  the  diffusion  term  and 
the  yield  at  low  fields  where  diffusion  is  already  most  significant. 
Crank  suggests  no  less  than  10  examples  of  concentration  dependent 
models  that  one  might  choose.  There  is  no  physical  reason  for  pre¬ 
ferring  one  of  them  over  another,  but  they  all  would  produce  a  cor¬ 
rection  in  the  right  direction  here. 


In  figure  32,  the  results  of  the  calculation  are  compared  with  the 
measurements  for  o  particles  incident  at  6  deg.  In  the  exposures  at  6 
deg,  the  normal  component  of  the  field  is  only  about  1/10  of  the  paral¬ 
lel  component.  Por  this  reason,  the  cylinders  move  apart  more  slowly 
than  they  move  past  each  other,  and  there  is  a  high  recombination  rate 
for  a  relatively  long  time,  at  least  where  the  charge  distributions 
still  overlap.  In  this  case,  the  effects  at  the  end  of  the  column  of 
charge  have  to  be  specifically  accounted  for  because  most  of  the  posi¬ 
tive  charge  that  escapes  recombination  is  near  the  Si-Si02  interface. 

N(t)  for  a  total  field  of  0.94  MV/cm  is  plotted  on  a  log  t  scale 
(fig.  33)  and  on  a  linear  t  scale  (fig.  34).  The  time  scale  can  be 
converted  to  a  ixjsition  scale  if  one  multiplies  by  The 
charge  near  the  Si-Si02  interface  contributes  more  to  the  observed 
threshold  voltage  shift.  This  effect  can  easily  be  accounted  for,  and 
the  experimental  results  in  figure  32  have  been  corrected  accordingly. 
Even  so,  one  can  see  that  the  measured  yield  is  still  higher  than  the 
calculated  yield.  The  probable  explanation  is  that  some  form  of  concen¬ 
tration  dependent  diffusion  is  again  appropriate. 


The  results  of  the  proton  calculation  are  compared  with  the  results 
of  the  experiment  in  figure  35  for  an  assumed  column  radius  of  3.5  nm. 
Once  again,  the  agreement  is  reasonably  good.  The  results  of  all  the 
calculations  are  summarized  in  table  5,  in  which  T  >  (L/cos  6)/ii_E  and 
the  yield  is  ( 1/NgT)/g  N(t)  dt. 
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Figure  15.  Initial  (time  =  0)  distribution  of  both 
positive  and  negative  charges  in  finite 
difference  solution  to  "exact"  Jaff^ 
equation  for  perpendicular  field  =  10^ 

V/cm. 


Figure  16.  Positive  and  negative  charge  densities  at 
time  *  10”^®  s  in  finite  difference  calcu¬ 
lation  for  perpendicular  field  =10®  V/cm. 


Figure  17.  Positive  and  negative  charge  densities  at 
time  =  10“^^  s  in  finite  difference  calcu¬ 
lation  for  perpendicular  field  =  10^  V/cm. 


DISTANCE  (mn) 

Figure  18.  Positive  and  negative  charge  densities  at 
time  =  10~15  s  in  finite  difference  calcu¬ 
lation  for  perpendicular  field  •  10®  V/cm 
(fig.  17  replotted  on  expanded  scale). 
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Positive  and  negative  charge  densities  at 
time  =  10“^^  s  for  perpendicular  field  = 
10®  V/cm  In  finite  difference  calculation* 


DISTANCE  (M*) 


Figure  22.  Positive  and  negative  charge  densities  at 
time  ^  3  X  s  for  perpendicular  field 

-  10**  V/cm  in  finite  difference  calcula¬ 
tion  (shows  importance  of  diffusion  at 
low  fields). 


DISTANCE  (Rin) 

Figure  23 •  Positive  and  negative  charge  densities  at 
time  =  10“^ 3  s  for  perpendicular  field  = 
105  V/cm. 
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Figure  24.  Positive  and  negative  charge  distributions 
at  time  ■  s  for  perpendicular  field 

■  2  X  10^  V/cm. 
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Figure  25.  Positive  and  negative  charge  densities  at 
time  =  10~^3  s  for  perpendicular  field  =  3 
X  10^  V/cm. 
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Figure  26.  Positive  and  negative  charge  distributions 
at  time  =*  10*’^^  s  for  perpendicular  field 
=  5  X  10^  V/cm. 
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Figure  27.  Positive  and  negative  charge  distributions 
at  time  =  s  for  perpendicular  field 

=  7  X  10®  V/cm. 


Figure  28.  Positive  and  negative  charge  distributions 
at  time  =  3  x  10”^**  s  for  perpendicular 
field  =  1.5  X  10®  V/cm. 
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Figure  29.  Positive  and  negative  charge  distributions 
at  time  =  3  x  10“l‘*s  for  perpendicular 
field  =  2  X  10^  V/cm. 


Figure  30.  Positive  and  negative  charge  distributions 
at  time  =  3  x  10”^**  s  for  perpendicular 
field  =  3  X  10^  V/cm. 


Figure  31*  Experimental  results  for  a 
particles  at  45  deg  (dots) 
repeated  along  with  results  of 
finite  difference  calculation 
(solid  line)* 


Figure  32.  Experimental  results  for  a 
particles  at  6  deg  (dots) 
compared  with  results  of 
finite  difference  calculation 
(solid  line)* 
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Figure  33.  Fraction  of  charge  surviving  at  time  t  for 
perpendicular  field  =  10®  V/cm. 


Figure  34.  fraction  of  charge  surviving  at  time  t 
for  perpendicular  field  “  10®  V/cm  (fig. 
33  replotted  on  linear  scale). 
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Figure  35. 

Comparison  of  proton  experi¬ 
mental  results  with  results  of 
finite  difference  calculation 
for  column  radius,  b,  =3  and 
3.5  nm. 
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TABtE  S. 

SUMMARY 

or  ETMITE  DIFFERESICE 

CALCOLATIOS 

RESULTS 

Unit 

Normal 

field 

(MV/cm) 

Total 

field 

(MV/cm) 

Time 

(s) 

r 

1  N(t)  dt 

Jo 

Fractional 

yield 

a  particle, 

0.1 

0.14 

4.73  K  10**^ 

9.17  X  10*^ 

0*0146 

45  deg 

0.2 

0.28 

2.36  «  10"'^ 

7.71  X  10** 

0.0246 

0.3 

0.42 

1.58  X  10"'2 

7.01  X  10** 

0.0334 

0.5 

0.71 

9.45  K  ID”*’ 

6.14  X  10** 

0.0489 

0.7 

0.99 

6.75  X  10"*3 

5.59  X  10** 

0.0623 

1.0 

1.41 

4.73  X  10'>3 

5.06  X  10** 

0.0804 

1.5 

2.12 

3.15  X  10"‘5 

4.42  X  10** 

0,1055 

2.0 

2.83 

2.36  X  to**’ 

4.01  X  10** 

0. 1278 

3.0 

4.24 

1.58  X  10"‘3 

3.49  X  10** 

0.1661 

o  particle, 

0.1 

0.94 

3.59  X  10”*5 

0.996  X  10** 

0.0208 

6  deg 

0.2 

1.89 

1.79  X  10"*3 

0.822  X  10** 

0.0345 

0.3 

2.83 

1.19  X  10"'^ 

0.729  X  10** 

0.0461 

Proton, 

0.1 

0.14 

4.73  X  10"'* 

6.43  X  10** 

0.0458 

45  deg 

0.2 

0.28 

2.36  X  10">2 

5.21  X  10** 

0.0743 

0.3 

0.42 

1.58  X  10"'2 

4.60  X  10** 

0.0981 

0.5 

0.71 

9.45  X  10"'5 

3.8?  X  10** 

Q.1378 

0.7 

0.99 

6.75  X  10*'’ 

3.42  X  10** 

0.1707 

1.0 

1.41 

4.73  X  10**’ 

2.98  X  10** 

0*2125 

1.5 

2.12 

3.15  X  10**’ 

2.50  X  10** 

0.2677 

2.0 

2.83 

2.36  X  10*'’ 

2.19  X  10** 

0.3130 

3.0 

4.24 

1.58  X  10*'’ 

1.80  X  10** 

0.3830 

5.  DISCUSSION 


This  work  shows  that  the  Jaff4  model  can  be  used  to  calculate  the 
ionization  of  Si02  by  heavy  charged  particles.  The  Jaffd  model  has  two 
"free"  pareuneters,  Ng  and  b,  but  Ng  can  be  determined  independently. 
This  work,  in  effect,  fixes  b  for  a  particles  and  protons  in  Si02  at 
~3.5  nm.  In  figure  35,  the  proton  experimental  data  are  plotted  along 
with  the  calculations  for  b  =  3  and  3.5  nm  to  give  an  idea  of  the  sensi¬ 
tivity  of  the  calculation  to  variation  in  b.  Similarly,  in  figure  36, 
the  45-deg  a  particle  data  are  replotted  along  with  calculated  yield 
curves  for  b  =  3  and  4  nm  as  well  as  b  «  3.5  nm.  One  can  see  that  the  b 
»  3.5  t  0.5  nm  is  a  reasonable  estimate  of  the  error  in  the  measurement 
of  b. 


One  reasonably  might  ask  why  b  should  be  in  the  range  of  3  to  4  nm; 
a  significant  body  of  literature  bears  on  the  gues> ion.  When  a  charged 
particle  passes  through  a  medium,  the  dominant  energy  loss  mechanism 
is  the  production  of  plasmons,  which  subsequently  decay  to  electron- 
hole  pairs  .**  ® '**  ®  Fbr  SIO2,  the  plasmon  energy  is  approximately 


Figure  36.  Comparison  of  a  particle  ex¬ 
periments  with  calculations 
for  column  radius,  b,  *  3, 
3.5,  and  4  nm  (indicates  how 
sensitive  results  are  to  vari¬ 
ation  in  b) . 


22  eV  (o)^  =  4iine2/m),  and  the  electron-hole  pair  energy  is  approximately 
18  eV.‘^^  This  4  eV,  which  Is  lost  In  thermallzatlon,  apparently  Is  lost 
through  the  emission  of  optical  phonons  having  an  energy  on  the  order  of 
0.1  eV.®®”53  The  mean  free  path  for  phonon  emission  by  hot  electrons  Is 
on  the  order  of  0.1  nm.^**  In  other  words,  the  total  path  length  Is  on 
the  order  of  4  nm  (40  events  that  are  0.1  nm  apart).  But  since  the 
electron-hole  pair  Is  on  a  random  walk,  the  effective  thermallzatlon 
radius  Is  on  the  order  of  times  0.1  nm  or  perhaps  0.6  nm,  clearly 
less  than  3.5  nm.  However,  If  one  takes  the  uncertainty  In  the  thermal 
energy  of  the  charges  to  be  on  the  order  of  kT,  one  can  calculate 
(Ap)^/2m  ~  kT  or  Ap  ~  4  X  10“^^  g-cm/s.  Then  from  the  uncertainty 
principle,  ApAx  >  h  or  Ax  >  2.4  nm.  Of  course,  the  uncertainty 
principle  gives  only  an  order  of  magnitude  estimate,  so  one  should 
probably  conclude  from  the  argument  only  that  b  will  fall  somewhere 
about  10”^  or  10“^  cm.  Certainly,  the  experimentally  determined  value  b 
=  3.5  nm  Is  within  this  range. 

Similarly,  one  can  estimate  the  spatial  extent  of  a  plasmon.  A 
plasmon  Is  basically  a  k  =  0  excitation,  except  that  the  quantum 
mechanical  uncertainty  Is  k^j_n  “  4nne^  =  m£(22  eV).®^  One  can  then 
calculate  AXj^^j^  ~  3  x  10”®  cm  or  roughly  an  order  of  magnitude  less  than 
the  extent  of  the  electron  wave  function.  Thus,  the  extent  of  the 
cylinders  of  charge  seems  to  be  determined  by  the  extent  of  the  electron 
wave  functions.  In  other  words,  the  radius  of  the  cylinders  Is  zero 
except  that,  quantum  mechanically,  nothing  Is  zero  to  arbitrary 
accuracy. 

Actually,  It  Is  somewhat  surprising  that  the  analysis  based  on  equa¬ 
tion  (15)  works  out  as  well  as  It  does  since  the  problem  should  really 
be  treated  quantum  mechanically.  Equation  (15)  Is  descended  directly 
from  the  Boltzmann  transport  equation,  which  Is  usually  described  as 
semlclasslcal  and  which  assumes,  for  example,  that  the  duration  of  the 
collisions  Is  small  compared  with  the  mean  free  time  between  colli¬ 
sions.  If  the  field  Is  taken  to  be  10®  V/cm,  then  the  drift  velocity 
saturates  at  ~2  x  10^  cm/s.  If  the  mean  free  path  Is  taken  to  be  0.1  nm 
again,  then  the  mean  free  time  Is  on  the  order  of  10~^®  s.  It  Is 
difficult  to  see  how  the  scattering  time  could  be  negligibly  small 
compared  with  10”^®  s. 

In  addition,  thermallzatlon  Is  assumed  to  be  Instantaneous,  but  the 
actual  thermallzatlon  time  Is  probably  on  the  order  of  10”^®  s.  For 
times  less  than  the  thermallzatlon  time  (that  Is,  for  hot  electrons),  D 
Is  probably  larger  than  the  value  given  below,  and  )i  Is  somewhat  smaller 
because  of  carrler-to-carrler  collisions. 

In  addition,  plasma  effects  such  as  the  shielding  of  external  fields 
are  neglected,  although  at  very  high  densities  (~10^®  cm”®)  these 
effects  should  be  significant. 
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In  spite  of  these  shortcomings,  the  theory  seems  to  work  reasonably 
well  at  predicting  the  yield  of  charge  In  the  limit  as  t  becomes 
large.  The  main  reason  seems  to  be  that  the  recombination  proceeds  very 
rapidly  at  early  times.  Most  of  the  charge  has  already  recoiid>lned 
before  the  electrons  thermalize  and  before  the  external  field  has  had 
much  effect.  The  screening  effects  become  less  significant  as  the  den¬ 
sity  decreases  (that  is,  the  assumptions  of  the  theory  become  better 
very  rapidly  as  time  Increases).  One  should  also  note  that  the  final 
yield  does  not  depend  on  the  precise  value  of  u  since  both  0  and  a  are 
assumed  proportional  to  p. 


The  classical  theory  used  here  consistently  underestimates  the  yield 
when  the  normal  component  of  the  field  is  small,  both  in  the  45-deg 
experiments  and  at  6  deg  and  for  both  a  particles  and  protons.  A  more 
complete  theoretical  treatment  with  a  larger  effective  D  to  account  for 
hot  electron  effects  would  probably  Improve  the  agreement  with  experi¬ 
ment. 


A  certain  amount  of  work  has  been  done  to  develop  a  full  quantum 
mechanical  treatment  of  transport  in  semiconductors and  in  dielec- 
trics.^^ However,  the  semlclasslcal  model  that  we  use  here  is 
empirically  effective  in  explaining  the  experimental  results,  even 
though  we  are  straining  the  assumptions  on  which  the  model  rests.  It  is 
interesting  how  often  a  theoretically  shaky  semlclasslcal  model  gives 
physically  reasonable  results  to  quantum  mechanical  problems — the 
Rutherford  scattering  cross  section  for  electrons  and  the  Bohr  atomic 
theory  are  two  other  examples. 


The  experiments  and  the  analysis  described  above  establish  that  the 
Jaffi  model  works  reasonably  well  for  o  particles  and  protons.  One  can 
easily  use  the  model  to  calculate  the  ionization  of  Si02  by  other  parti¬ 
cles  or  at  other  energies  since  these  factors  enter  the  calculation  only 
through  Nq,  the  number  of  electron-hole  pairs  per  unit  path  length.  In 
table  6,  we  present  calculated  results  for  the  cosmic  ray  nuclei  and 
^®Pe  along  with  our  typical  a  particle  and  proton  results.  The  energies 
chosen  for  and  ®”Pe  correspond  to  the  maximum  ionization  for  those 
particles,  and  the  a  particle  and  proton  energies  are  chosen  to  match 
our  experimental  conditions.  One  can  see  from  table  6  that  the  initial 
ionization  density  for  is  roughly  25  times  greater  than  for  an  a 
particle.  But  the  amount  of  charge  that  escapes  recombination  is  only 
about  twice  as  great  for  The  reason  is  that  the  recombination 
term,  -an^n-,  is  roughly  600  times  greater  for  at  first,  and  the 
extra  free  charge  is  eliminated  rapidly. 
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TABLE  6. 

CALCULATED  RESULTS 

FOR  PROTONS, 

a  PARTICLES,  12c, 

AND 

Particle 

Particle 
energy,  B 
(MeV) 

Energy  loss 
dB/dx 

( MeV/g/cra^ ) 

,  Electron 

hole 

(pairs/cro) 

Fractional  yield 
(field  »  1.4  MV/cm) 

Yield 

(holes/cm) 

‘h 

0.7 

243 

2.97  X  10^ 

0.2120 

6r30 

X  10® 

•♦He 

2.0 

1,068 

1.33  X  10® 

0.0885 

1.18 

X  10^ 

12c 

4.0 

5,270 

6.44  X  10® 

0.0242 

1.56 

X 

o 

100 

30,680 

3.75  X  10® 

0.0060 

2.25 

X  10^ 

6.  IMPLICATIONS  FOR  MICROELECTRONIC  DEVICES 

The  possibility  of  permanent  failures  of  electronic  devices  from 
single  particles  is  a  timely  question  in  view  of  the  current  interest  in 
temporary  upsets. One  might  reasonably  ask  how  small  the  devices 
would  have  to  be  before  permanent  errors  would  be  observed.  Recently, 
Srour  et  al®^  considered  the  possibility  that  a  single  neutron  might 
cause  enough  displacement  damage  to  ruin  a  submicrometer  device.  In 
addition,  Oldham  and  McGarrity®^  have  presented  a  worst  case  calculation 
indicating  that  a  single  a  particle  might  cause  enough  ionization  to 
cause  a  device  to  fail. 

This  worst  case  analysis  rests  on  the  assumptions  that  ( 1 )  recombin¬ 
ation  can  be  neglected  and  (2)  all  the  charge  reaching  the  Interface  is 
trapped  there.  The  assumption  of  100-percent  trapping  is  reasonable  in 
some  cases, but  the  assumption  of  100-percent  yield  is  clearly 
unjustified  in  view  of  the  experimental  results  presented  above.  For 
100-percent  yield,  the  device  discussed  by  Dennard  et  al®**  in  their 
well-known  scaling  law  paper  ( 1  x  1  ym  with  35-nm-thick  Si02)  would 
exhibit  a  threshold  voltage  shift  (&Vq>)  of  about  100  mV  after  being 
struck  by  a  2-MeV  a  particle.  However,  if  recombination  were  taken  into 
account,  the  yield  would  be  only  6  to  10  percent  in  the  field  range  of 
1  to  2  MV/cm,  and  AVij>  would  be  only  about  6  to  10  mV.  Many  present  day 
circuits  fail  when  AVq.  =  100  mV,  but  they  can  withstand  AVip  =10  mV. 

One  can  calculate  AV>p  as  a  function  of  device  dimensions  starting 
from  equation  (6).  If  the  total  charge  produced  by  an  ionizing  par¬ 
ticle  is  given  by  equation  (6),  then  the  charge  escaping  recombination 
is  Q  =  f(E)Qg,  as  before,  and  the  charge  trapped  at  the  interface  is 
f^f(E)Qg,  where  fip  is  the  fraction  of  holes  that  reach  the  Si-Si02 
Interface  and  are  trapped  there.  Then 

'"'t  ■  ^  ^  ^  “HP* 


(28) 


where  A  Is  the  active  area  of  the  device*  If  we  assume  a  2-MeV  a 
particle  with  f{E)  =  10  percent  (that  is,  E  *  2  MV/cm)  and  fip  =  100 
percent,  we  can  plot  equation  (28)  for  devices  of  different  areas  as  a 
function  of  SIO2  thickness  with  the  results  shown  In  figure  37.  The 
results  in  this  plot  do  not  Include  short  channel  radiation  effects,®® 
which  are  poorly  understood  but  which  may  be  important.  In  considering 
figure  37,  one  should  remember  that  voltages  are  scaled  down  somewhat  as 
devices  are  reduced  in  size.  Fbr  this  reason,  tolerable  also  will 
tend  to  shrink  in  the  future.  On  the  other  hand,  we  have  assumed  ffj>  = 
100  percent,  which  is  realistic  only  for  unhardened  commercial  8102*  Eor 
hardened  Si02^  the  trapping  fraction  can  be  reduced  by  a  factor  of  10  or 
perhaps  more,  and  AV^  can  be  reduced  by  the  same  factor. 

Actually,  the  calculation  of  AVi^  in  figure  37  is  misleading  because 
we  are  dividing  total  trapped  charge  by  total  capacitance.  In  effect, 
we  are  assuming  that  the  charge  will  be  trapped  uniformly  across  the 
active  area  of  the  device,  but  one  would  really  expect  a  nonuniform 
spatial  distribution  of  charge  as  indicated  schematically  in  figure  38. 


Figure  37.  Calculated  threshold  voltage 
shift  for  various  small  area 
devices  as  function  of  SIO2 
thickness.  (This  calculation 
assumes  2-MeV  a  particle  inci¬ 
dent  at  45  deg  emd  10-percent 
yield;  that  is,  applied  field  » 
2  MV/cm.  Fbr  lower  fields  or 
for  other  particles  [protons]  or 
other  energies,  threshold  volt¬ 
age  shift  is  smaller  than  shown 
here. ) 
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Figure  38.  Schematic  of  MOSFET  struck  by  a 
particle. 


One  can  do  a  fairly  simple  "particle  pushing"  calculation  to  esti¬ 
mate  the  spatial  distribution  of  the  charges  as  they  arrive  at  the 
interface.  We  distribute  randomly  the  correct  number  of  positive 
charges  in  approximately  the  correct  volume  and  calculate  the  coulomb 
repulsion  between  each  pair  of  charges.  Then  following  the  hole 
transport  model  of  McLean, we  let  each  hole  hop  1  nm  parallel  to  the 
total  field  that  it  sees.  Then  we  stop  the  particles,  recalculate  the 
fields,  and  let  the  charges  hop  again.  This  process  is  continued  until 
all  the  charges  are  trapped  at  the  interface.  The  charges  generated 
near  the  interface  reach  it  first  and  are  trapped  near  the  exit  point  of 
the  ionizing  particle  from  the  Si02.  As  these  charges  build  up,  they 
generate  a  field  that  tends  to  cancel  the  applied  field.  Charges 
arriving  later  have  to  spread  from  the  particle  track  until  they  are 
outside  the  central  space  charge  region,  and  only  then  do  they  move  to 
the  interface  under  the  influence  of  the  applied  bias. 

This  process  is  indicated  schematically  in  figure  39  for  an 
obliquely  incident  particle.  The  innermost  zone  "fills  up"  first  and, 
as  time  passes,  the  remaining  zones  fill  up  one  at  a  time  in  order  of 
increasing  distance  from  where  the  particle  hit,  r. 

Figure  40  shows  the  results  for  an  a  particle  normally  incident 
on  SIO2  100  nm  thick,  where  the  applied  field  is  10^  V/cm  and  recom¬ 
bination  is  neglected.  Even  though  the  density  at  r  <  10  nm  is  very 
high,  the  total  charge  in  the  first  zone  is  only  53  holes  in  this 
example.  These  53  holes  represent  a  density  greater  than  10^^ 
charges/cm^ ,  corresponding  to  an  immense  field.  These  holes  were 
created  very  near  the  interface,  and  the  field  that  they  generate 
prevents  any  more  charges  from  transporting  into  the  r  <  10  nm  zone. 
Charges  arriving  later  at  the  Interface  have  to  move  out  to  large  r 
before  the  applied  field  can  drive  them  to  the  interface. 
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Figure  39.  Qualitative  illustration  of  path  followed 
by  holes  hopping  to  interface  in  particle 
pushing  calculation. 
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Figure  40.  Charge  density  at  interface, 
assuming  normally  incident 
2-MeV  a  particle.  lOO-nm- 
thick  Si02«  applied  field  > 
10®  V/cm,  and  100-percent 
yield.  (Dashed  line  shows 
average  [mean]  charge  dens- 


The  charge  distribution  that  results  from  a  calculation  of  this  kind 
depends  greatly  on  the  yield*  In  figures  41  and  42,  results  are  given 
for  Aqx  =  20  nm,  where  the  field  is  10^  V/cm.  The  a  particle  is 
incident  above  the  origin  and  exits  at  x  =  20  nm  and  y  =  0.  The  figures 
are  what  one  would  see  looking  down  from  above  on  the  interface,  where 
each  dot  represents  a  trapped  hole.  In  figure  41  the  yield  is  100 
percent,  but  in  figure  42  it  is  only  10  percent.  The  maximum  density  is 
much  greater  for  the  higher  yield,  and  the  size  of  the  "footprint"  also 
is  much  greater. 

A  more  realistic  case  for  electronic  devices  today  and  in  the  near 
future  is  shown  in  figure  43.  The  splatter  pattern  in  figure  43  is  for 
an  a  particle  incident  at  45  deg  on  50-nm  Si02,  where  E  =  10^  V/cm  and 
the  yield  is  assumed  to  be  10  percent.  As  before,  the  particle  is 
incident  above  the  origin,  and  the  exit  point  is  x  =  50  nm  and  y  =  0. 
One  can  see  that  the  footprint  is  on  the  order  of  100  nm  in  diameter, 
which  is  much  less  than  the  size  of  a  1-pm  device  or  even  a  0.5— pm 
device.  If  one  imagines  this  charge  distribution  placed  in  a  0.5-iim 
channel,  one  can  calculate  a  local  flatband  voltage  shift  from  the  local 
charge  density.  There  would  be  a  large  local  shift,  but  most  of  the 
channel  would  not  be  affected  at  all.  This  result  is  illustrated  in 
figure  44. 


■  (■Ml) 

Figure  41.  Charge  density  at  interface,  assuming 
that  2-MeV  a  particle  struck  20-nm- 
thick  Si02  at  45  deg  with  applied 
field  >•  10^  V/cm  and  100-percent 

yield.  (Bach  dot  represents  one 
trapped  hole.  Figure  shows  distri¬ 
bution  of  holes  that  one  would  see 
looking  down  from  above  Si-Si02 
interface. ) 
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Figure  42.  Charge  density  at  Interface,  assuming 
that  2-MeV  a  particle  struck  20-nm- 
thick  SiOj  at  45  deg  with  applied 
field  =  10^  V/cm  and  10-percent 

yield.  (Each  dot  represents  one 
trapped  hole.  Figure  shows  how  spot 
size  depends  on  yield. ) 
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Figure  43. 

Charge  density  at  interface, 
assuming  that  2-MeV  a  particle 
struck  50-nm-thick  Si02  at  45 
deg  with  applied  field  =  10^ 
V/cm  and  10-percent  yield. 
(Each  dot  represents  one 
trapped  hole.  Si02  thickness 
is  more  nearly  typical  of 
present  day  devices. ) 
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Figure  44.  Local  flatband  voltage  shift  caused 
by  placing  charge  distribution  in 
figure  43  in  middle  of  device  with 
0.5-pm  channel  length. 


Thus,  it  seems  unlikely  that  a  single  charged  particle  will  cause  a 
threshold  voltage  shift  large  enough  to  cause  a  total  dose  failure  in  a 
device,  at  least  until  devices  are  scaled  down  well  below  1  um.  How¬ 
ever,  the  small  charge  blobs  such  as  in  figure  43  could  cause  unforeseen 
reliability  or  stability  problems  such  as  punch  through  or  Increased  hot 
electron  injection.  One  should  remember  that  highly  integrated  circuits 
are  increasingly  difficult  to  design  and  build  even  in  the  absence  of 
radiation.  Adding  a  few  unexpected  blobs  of  charge  to  these  circuits  is 
a  little  bit  like  throwing  a  handful  of  sand  into  a  jet  engine.  The 
effects  are  difficult  to  predict,  but  they  cannot  possibly  be  benefi¬ 
cial. 


7.  SUMMARY  AND  CONCLUSIONS 

We  have  measured  the  recombination  of  charge  produced  in  Si02  by  two 
kinds  of  heavy  charged  particles,  a  particles  and  protons.  The  fraction 
of  charge  that  escapes  recombination  is  relatively  small  at  fields  on 
the  order  of  1  to  2  MV/ cm.  The  yield  for  a  particles  is  10  percent  or 
less,  and  the  yield  for  protons  is  25  percent  or  less.  On  the  other 
hand,  for  relativistic  electrons  the  yield  is  90  percent  or  more  at 
these  fields. 
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The  results  of  these  experiments  are  explained  reasonably  well  by 
the  Jaf?^  columnar  recombination  model  when  a  column  radius  of  3.5  nm  is 
assumed.  The  Jaff€  model  is  well  known,  but  the  column  radius 
for  Si02  had  not  previously  been  determined.  Since  the  column  radius  is 
presumably  a  property  of  the  material,  the  Jaff6  model  can  now  be  used 
to  calculate  the  yield  of  charge  for  other  particles  and  other  energies. 

Finally,  the  implications  of  the  charge  yield  measurements  for  the 
operation  of  microelectronic  devices  are  considered.  The  charge  that  is 
produced  by  an  a  particle  and  escapes  recotnbination  is  unlikely  to  cause 
a  total  dose  failure  in  a  device  until  the  devices  are  much  smaller  than 
1  pm.  However,  the  circuits  are  complicated  enough  that  the  charge  from 
a  single  a  particle  could  make  it  more  difficult  to  solve  some  relia¬ 
bility  problems. 
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The  FORTRAN  listing  of  the  finite  difference  code  is  given  below. 
Subroutine  INIT  initializes  certain  variables,  subroutine  INT  integrates 
the  total  charge  and  updates  the  yield  integral,  subroutine  LAPLAC  cal¬ 
culates  the  Laplaclan  operator,  subroutine  DNDX  calculates  the  deriva¬ 
tives  3n/3x,  subroutine  ONDT  calculates  the  derivatives  3n/3t,  subroutine 
NTKP1  calculates  the  new  charge  densities  for  the  next  cycle,  and  sub¬ 
routine  OOTPOT  writes  the  results.  The  main  progrcun  supervises  the  sub¬ 
routines  . 


IMPLICIT  R6AL»8  (A-H,0-2» 

COMMON  X(,  Y(IOI)  t  RHON(20l.lOI>»  RH3P  (  201 1 101 ) 

COMMON  □2KhQN(20lf 101) »  02RH0PI 201 . 101 ) ,  00XRMN(201 , 101 > 

COMMON  OOXRHP(201.10l) t  ONNOT (2C1 » ICl ) >  ONPOT ( 201 « 101 ) » TPRINT 1 10 ) 
COMMON  ciBtONfDPf AL?HA,RHOLfUN,UP,TOT.N*TOTP,DELTOtDeLT»  TRATIO 
COMMON  T,0ELX,0ELYtYIEL0tNPRINT,lPRlNT,K 

COMMON  NXMAXtNYMAXfNTMAXiNMAXtNXMXMltNYMXMly JMiNtJMAXtIPLAG 

CALL  INir 

CALL  INT 

CALL  OUTPUT 

DO  200  K=lfNTMAX 

CALL  LAPLAC 

CALL  ONOX 

CALL  ONDT 

CALL  NTKPl 

CALL  INT 

CALL  OUTPUT 

T»T40ELT 

OELT*TRATIO*OeLT 
IF  (IFLAG.GT.NPRINT)  GO  TO  201 

200  CONTINUE 

201  STOP 
END 
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SUBROUTINE  INIT 
implicit  REAL*8  JA-H,a-Z> 

COMMON  X(20ll.  YUOl)  t  RHONt  201.  ICl )  ,  RHOP  (201tlQl ) 

COMMON  02RH0N(201.101)*  02RHOP( 201 • 101 ) t  ODXRHN (2C1 » 101 1 
COMMON  00XRHP(2Cltl01) .  ONNOT (201 » 101 ) t  0NP0T(201»101) .TPRINT(IO) 
COMMON  EtBtONtOP«ALPHA,RHOLtUNyUP,TOTN*TOTPtOELTO»DELT»  TRATIO 
COMMON  T.DELXtOELYtYIELOtNPKINTtlPRINTyK 

COMMON  NXMAXfNYMAX.NTMAXtNMAXtNXMXMltNYMXMlf JMINtJMAXtIFLAG 

READ  (StlOO)  E.StON.DPtALPHA.RHOLtUNiUP 

READ  (StlOO)  OELXtOELYtOELTOtTRATIO 

READ  (St  101)  NXHAXtNYMAXtNTHAXtNPRlNTt JMlNt JMAX 

READ  (StlCO)  (TPRINT(I)tI»ltNPRINT) 

100  FORMAT  (9D10.3) 

101  FORMAT  (bllO) 

WRITE  (6tl02) 

102  format  (IHI) 

WRITE  (6tl03)  EtStONtOPtUNtUPtALPHAtRHOL 

103  FORMAT  (•  E»‘t010.3t5X,*  8=  ‘tOlO.StSXt*  DN=  •tD10,3t5Xt»  DP=  • 

X  010.3tSXt'  UNs  •,010.3t5Xt*  UP-  • t010.3tSXt//t *  ALPHA*  SOlO. 

X  SXt*  RHOL-  *,010.3) 

WRITE  (6,104)  DELXtOELYtOELTOtTRATlO 

104  FORMAT  (•  0ELX=*,010.3tSXt*  OELY*  *t010;3t5Xt*  OELTO*  'tOlO.StSX, 

X  •  TRAT10=*,010.3) 

WRITE  (6,105)  NXMAX,NYMAX,NTMAX,NPRINT,JMIN,JMAX 

105  FORMAT  (•  NXMAX=  •,I6,6X,*  NYMAX*  'tlGtbX,*  NTMAX=*,  16, 6X, 

X  •  NPR1NT=  •,16,6X,*  JMINs  •,16,6X,  *JMAX=  *,16) 

WRITE  (6,106)  (TPR1NT(I),1*1,NPRINT) 

106  format  (•  TPRINT*  *, 10(010.3, 2X ) ) 

NMAX«2'*r4XMAX-l 

DO  200  1*1,NMAX 
X(1)«(1-NXMAX)*DELX 

200  CONTINUE 
Y(l)*  0.00 

DO  203  J=2,NYMAX 
Y(J)=  (J-1)*0ELY 

203  CONTINUE 

00  201  J*1,NYMAX 
00  201  I*1,NMAX 
RSQ*-( X ( I ) **2+ Y( J (•♦Z 
IF  (RSO.LT.-ICO.OO)  GO  70-204 

RHONd,  J)*  (RH0L/(3.14159*B4B))*0EXP{-(X(I)**2*Y(J)**2)/B**2> 
RH0P(I,J)*RH0N(I,J) 

GO  TO  205 

204  RHON(I,J)*0.00 
RHOP( 1, J)*O.DO 

205  02RHON(I,J)*0.00 
02RH0P(1, J)*O.OC 
00XRHN(1,J)>0.00 
00XRHP(I,J)sC.D0 

201  CONTINUE 
1FLAG«0 
IPRINTsO 
TaO.OO 
OELT-DELTO 
YIELD  *0.00 
RETURN 

END 
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SUBROUTINE  LAPLAC 
IMPLICIT  HEAL*3  lA-HtO-ZI 

COMMON  X(201),  YdOll  »  RH0N(2Cltl01) »  RH0P(201»101) 
common  02KH0NI201tl01lt  nZRHOPi 201* 101 ) t  00XRHN<201 t 101 ) 

COMMON  00XRHP(201tl01)»  ONNOTI 201  *  101 ) «  DNPOT <201 1 101 ) > TPRINT ( 10) 
COMMON  E,B,ON,DP,ALPHA,RHOLtUN,UPtTOTN,70Ti»,OELTOtDELTt  TRATIO 
COMMON  TtOELXtOELYtYlELOfNPRINTtIPRINT.K 

COMMON  NXMAXfNYMAX«NTMAX»NMAX«NXHXMl*NYMXMlt JMlNtJMAXtlrLAG 

NXMXMl=NMAX-l 

NYMXMlsNYMAX-1 

02RH0N<1,1)=  (RH0N<2fl)«RH0N(l  »2 )-2*RH0N( 1 , 1 ) )/0ELX**2 
D2RH0P(l,l)a  (RH0P(2,l)>RHQP< 1  *2 )-2*RH0P (1 . 1 ) )/DELX**2 
DO  200  J=2tNYMXMl 

02RH0N(  It  J)  =  <RH0N(2*JdRH0N<l  *  J-1 )  ♦RHONC 1.  J^l  )-3*RH0N(  1 «  J)  )/ 

X  0£LX**2 

02RH0P< It Ji»(RH0P(2.J)«RH0Pll  t J-1) ♦RHOPC It J«^ll -3*RH0P( 1  * J ) )/ 

X  0ELX**2 

200  CONTINUE 

02RH0N ( 1 1 N YM AX ) s  <  RHON <  2  tN  YM AX ) ♦RHON 1 1  * N YMXM 1 » -2 *RHON ( 1 1 NYM AX )  ) / 

X  0ELX**2 

02KH0P< ltNYMAX»  =  (RHOP( 2 tNYMAX )>RH0PI1 1 NYMXM 1 ) -2*RM0P ( 1 1 NYMAX )  )/ 

X  0ELX**2 
00  201  1-2  tNXMXMl 

02RH0N(Itl)s(RH0N(I-»l,l)4RH0N|I-l  1 1)  «2*RH0N(  I  *2 ) -4*RH0N  ( 1 1 1 )  )/ 

X  0ELX**2 

D2RH0P(ltl)*<lM0P(I>ltU+RH0P<I-l  tl)^2*RH0PII,2»-4*RH0P(Itl))/ 

X  0ELX**2 

02RH3N(  I  tNYMAX  J*(RH0N(  I-»ltNYMAX)+RHON<  I-ltNYMAX»  ♦RHON  1 1  tNYMXMl  )- 
X  3*RH0P(I,NYMAX) )/JELX**2 

02RH0P  < I iNYMAX )  =  ( RHOP ( I ♦! tNYMAX )*RH0P ( 1-1 ,NYMAX ) ♦RHOP ( I tNYMXMl  )- 
X  3*RH0N(I,NYMAX) »/06LX**2 

201  CONTINUE 

DO  202  J=2tNYMXMl 
DO  202  I=2tNXMXMl 

02RH0N< I, JI=(RHON(I-l, J»>RhON<I^ltJ»*RHON(I t J+1 )+RHON( I t J-1 »- 
X  A*RHON(I,JI )/DELX*A2 

02RH0P(  It  JI=(RHaP(I-l,  J»  +  RHOP(  !♦!*  J»4RH0P(I  tJ^D^RHOPd  tJ-D- 
X  A*RHnP(Ifjn/0cLX**2 

202  CONTINUE 

0?RH0N(  NMAXt  l|s(RHaN<NXMXMltldRHON(  NMAX  t  2  I-2*RH0N  (  NMAXtl))/ 

X  0ELX**2 

02RHOP(  NMAXtllsIRHOPtNXMXMltD^RHrjPf  NMAX 1 2 )-2*RH0P( NMAX t  1))/ 

X  0ELX**2 

D2RH0N(  NMAX, NYM AX  I =< RHON (NXMXMl tNYMAX I^RHONI  NMAXtNYMXMl )-2*RH0N 
X  (  NMAXtNYMAX) )/0ELX*«2 

D2RH0P (  NMAX  tNYMAX »  =  <  RHOP (NXMXM 1 ,NYMAX  » ♦RHOP I  NMAX  tNYMXMl ) -2*RH0P 
X  (  NMAX,NYMAXn/0ELX**2 
DO  203  J*2 tNYMXMl 

DZRHONI  NMAXtJ)»  ( RHON  <  NXMXMl  ,  J  )^RHON  I  NMAX,  J-l)  ♦RHON  <  NMAXtJ^D- 
X  SPRMONI  NMAXtJ » J/0ELX**2 

02RH0P(  NMAXtJ)*  (RhOPINXMXM1,J)^RnOP(  NMAX , J-1 l♦RHOP(  NMAXtJ^ll- 
X  3*RH0P(NMAX,  J))/0ELX**2 

203  CONTINUE 
RETURN 
END 
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subroutine  ONOX 

IMPLICIT  REAL»6  (A-H,0-I» 

COMMON  X(201)t  Y(lOl)  t  RHONt 201 « 101 > t  RHOP (201 tlOl I 
COMMON  02RHON(201tl01)»  02RH0PI 20 1 , 10 1 1 1  00XRHN(2C1 *  101 ) 

COMMON  □i:)XRHP(201«101)  t  ONNDT  (  201  «10f  )  .  ONPOT ( 201 1 101 )  , TPRINT  ( 10 1 
COMMON  Et6*DN,0P,ALPHA,RH0LtUN,UP,T0TNtT0TPtDELT0fDELT,  TRATIO 
COMMON  T.DcLX.OELYtYIELO.NPRINT.IPKlNl.K 

COMMON  NXMAX,NYMAX,NTMAX*NMAX,NXMXMltNYMXMlt JMIN, JMAXtIFLAG 
DO  200  J=ltNYMAX 

DOXRHNd  t  J)s(RH0N(2tJl-RH0N(  It  J>  )/DELX 
DDXRHP(  1  ,J)  =  (RH0P(2,J»-RM0PI  1,  J)  )/DEH 
OOXRHNt  NMAXtJ)  =  (RHON(  NMAX.  J  )-KHC’N  (NXMXMl  *  J  )  l/OELX 
00XRHP(  NMAXt J)=(RNOP(  NMAX t J l-RHOP (NXMXMl , J M/OELX 
200  CONTINUE 

00  202  J=ltNYMAX 
DO  201  1S2.NXMXM1 

OOXRHN ( 1 , J )  =  ( RmON ( I ♦ 1 1  J 1-RHON ( 1 -1 1  J  t ) / ( 2*DE  LX ) 

DDXRHP( I,J)=(RM0P(I*1, J)-SHOP(1-1,J»)/(2*OELX) 

2C1  CONTINUE 
202  CONTINUE 
RETURN 
END 


SUBROUTINE  DNOT 
IMPLICIT  REAL*3  (A-H,0-Z1 

COMMON  X(201)t  Y(lOl)  ,  RHONI 201t 101 1 t  RHOP (201 « lOl 1 
COMMON  02RM0N(20Itl01) t  02RH0P( 201t 101 ) t  DOXRHN(201 t 101 ) 

COMMON  0OXRMP(20ltlont  ONNDT  ( 201 1 101  >  »  DNPOT  (201t  101 )  t  TPRINT  ( 10 » 
COMMON  EtBtONtOP, ALPHA, RHOL»UN,UPtTOTNtTOTP,DELTOtOELTt  TRATIO 
COMMON  TtCELXtOELY, YIELD, NPRINTtIPRINTtX 

COMMON  NXHAX,NYMAX,NTMAX,NMAX,NXMXH1,NYMXM1, JMIN,JMAX,IFLAG 
00  201  J=l,NYMAX 
00  200  I=1,NMAX 

ONNOT( I , J|=0N*02RH0N(I ,J»^UN*E*DDXRMN(1, JI-ALPHA*RHON(I,JI* 

X  RHOP(I,J) 

ONPOT(  I ,J»=«0P*02RH0P( I ,J)>UP*E*ODXRHP( I,J)-ALPHA*RH0N(1,J)P 
X  RHOPdtJI 

200  CONTINUE 

201  CONTINUE 
RETURN 
END 
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subroutine  NTRPl 
IMPLICIT  REAL*8  (A>H,0-Z) 

COMMON  XIZOlt.  YClOl)  t  RHONlZCltlOl) «  RHOP 1201 , 101 ) 

COMMON  02RH0NI201il01)t  02RH0PI 201 « 101 ) «  D0XRHN(201, ICl I 
COMMON  0DXRhP(2Cl,101) ,  0NNOT«201,101 ) ,  ONPOT (201 1 101 ) » TPRINTC 10 1 
COMMON  Etb.ON.OP.ALPnA,RHDLtUN,UP«rOTN,TOrP«OELTOfOELTt  TRATIO 
COMMON  TtOELX.OELYt YIELD, NPRINT,I*PRINT,K 

COMMON  NXMAX,NVMAX,NTMAX,NMAX,NXMXM1,NYMXM1, JMlNtJMAXtlFLAG 
00  201  J«l,NYMAX 
00  200  1*1,  NMAX 

RHONC 1 , J I =RHON ( 1 , J  J ♦ONNOT ( I , J 1*0ELT 
RHOP(  1,J)  sRhOPd.J)  ♦ONPOT  II,  JI*0ELT 

200  CONTINUE 

201  CONTINUE 
RETURN 
END 


subroutine  INT 
IMPLICIT  RcAL*6  <A-H,0-ZI 

COMMON  X(20U,  YllOll  ,  RH0NI2C1, 101)  ,  RHO®  1201 , 101 ) 

COMMON  02RH0N(201,10U,  02RHOP(  201 , 101  >  ,  D0XRHN(201 , 101 ) 

COMMON  oaXRHP(20l,igil ,  ONNOT <201,101 ) ,  ONPOT<201,101),TPR1NTIIO> 
COMMON  s,b,dn,op,alpma,rhol,un,up,totn,totp,oelto,delt,  TRATIO 

COMMON  T,0£LX,D=LY,YI£L0,NPRINT,IPR1NT,K 

COMMON  NXMAX, NVMAX,NTMAX, NMAX, NXMXM1,NYMXM1,JMIN, UMAX, IFLAC 

TOTNsO.OO 

TOTPsO.DO 

00  200  1=1, NMAX 

TOTN»  RH0NH,1)»DELX*0ELY^T0TN 
TOTP=  RHOPd,  l)>*DeLX*OELY+TOTP 

200  CONTINUE 

00  201  J=2,NYMAX 
00  201  1=1, NMAX 

T0TN»2*RhQN( I , J)*DELX*0ELY+T0TN 
T0TP=2*RH0Pd,  J)*0ELX*0ELY>T0TP 

201  CONTINUE 
YIEL0=YIEL0+T0TP»DELT 
RETURN 

END 
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SUBROUTINE  OUTPUT 
IMPLICIT  REAL*8  (A-M,0-Z» 

COMMON  X(20I).  Y(lOl)  t  RHON( 201. 101 > »  RHOP (201 . 101 ) 

COMMON  D2RH0N(201tl01) >  02RH0P( 201 . 101 ) «  00XRHN(201> 101 ) 

COMMON  00XRHP(2C1.101lt  ONNOT(2Cltl01 ) .  ONPOT (201 . 101 ) . TPRINT ( 10) 
COMMON  E.B.ON.OP.ALPHA.RHOL.UN.UP.TOTN.TOTP.DELTO.DELT,  TRATIO 
COMMON  T.DELX.OELY.YIELO.NPRINT.IPRINT.K 

COMMON  NXMAX.NYMAX.NTMAX.NMAX.NXMXMl.NYMXMl , JMlN.JMAXtlFLAG 
WRITE  (6.101)  K.T.TOTN.TOTP. YIELD 

101  FORMAT  (•  CYCLE  NO*  *15. 5X.*  T*  '.OIO.O.SX.*  TOTN»  *.010.3. 

X  5X.*  TOTPs  ‘.OlO.a.SX.*  YIELD*  *.010.3) 

IF  (T.GE.TPRINT(IFLAG)  )  1PR1NT=1 
IF  (IPRlNT.EQ.l)  GO  TO  201 
GO  TO  202 

201  WRITE  (6.102) 

102  FORMAT  (IHI) 

00  203  1*1.NMAX 

WRITE  (6.103)  ( (X( 1 ).Y(J).RH0N(1.J).RH0P(I*J)  ).J«JMIN.JMAX) 

103  FORMAT  (•  X=* .010.3. 5X.*  Y** . 01 0.3. 5X . •  N=  *.010.3.SX. 

X  •  P=  *.010.3) 

203  CONTINUE 
IPRINT=0 
1FLAG=IFLAG*1 

202  RETURN 
END 
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